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Linearly  forced  isotropic  turbulence 

By  T.  S.  Lundgren  f 


Stationary  isotropic  turbulence  is  often  studied  numerically  by  adding  a forcing  term 
to  the  Navier-Stokes  equation.  This  is  usually  done  for  the  purpose  of  achieving  higher 
Reynolds  number  and  longer  statistics  than  is  possible  for  isotropic  decaying  turbulence. 
It  is  generally  accepted  that  forcing  the  Navier-Stokes  equation  at  low  wave  number  does 
not  influence  the  small  scale  statistics  of  the  flow  provided  that  there  is  wide  separation 
between  the  largest  and  smallest  scales.  It  will  be  shown,  however,  that  the  spectral 
width  of  the  forcing  has  a noticeable  effect  on  inertial  range  statistics.  A case  will  be 
made  here  for  using  a broader  form  of  forcing  in  order  to  compare  computed  isotropic 
stationary  turbulence  with  (deca3ung)  grid  turbulence.  It  is  shown  that  using  a forcing 
function  which  is  directly  proportional  to  the  velocity  has  physical  meaning  and  gives 
results  which  are  closer  to  both  homogeneous  and  non-homogeneous  turbulence. 

Section  1 presents  a four  part  series  of  motivations  for  linear  forcing.  Section  2 puts 
linear  forcing  to  a numerical  test  with  a pseudospectral  computation. 


1.  Motivation  for  Linear  Forcing 

1.1.  Linearity  of  energy  production  in  non-homogeneous  turbulence 
In  shearflow  turbulence  the  equation  for  the  fluctuating  part  of  the  velocity,  u',  is 



-f  u • Vu'  + u'  • Vu  -I-  u'  • Vu'  — V • u'u'  = -Vp'lp  -t-  z/V^u'.  (1) 

at 

The  third  term  on  the  left,  u'  • Vu  appears  in  the  turbulent  energy  equation  as  the  energy 
production  term,  < u'  • Vu  • u'  >.  (Both  angle  brackets  and  over  bars  are  used  to  denote 
averages.)  In  (1)  it  appears  as  a forcing  term  proportional  to  u'.  This  suggests  that  for 
isotropic  homogeneous  turbulence  it  might  be  appropriate  to  force  a stationary  flow  with 
a driving  term  proportional  to  the  velocity.  Of  course,  for  isotropic  turbulence  there  is 
no  mean  velocity  gradient,  so  the  only  way  to  have  the  flow  be  isotropic  and  stationary 
is  to  have  the  forcing  be  isotropic.  It  is  proposed  here  to  use 

^ + u • Vu  = — Vp/p-l- i/V^u-l- f (2) 


with  the  driving  force 


f = Qu, 


(3) 


where  Q is  a constant.  The  prime  on  the  velocity  has  been  omitted  here  and  henceforth 
with  the  understanding  that  < u >=  0.  The  turbulent  energy  equation  is  now 


1 5 < u • u > 

2 dt 


-e-i-Q  < u-u  >, 


(4) 


where 


e ==  — z/  < u • V^u  > 


f University  of  Minnesota 


(5) 
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is  the  mean  dissipation  rate  and  the  last  term  could  be  called  isotropic  turbulent  produc- 
tion. For  stationary  turbulence  therefore  (setting  the  time  derivative  to  zero) 

Q = e/3ff2  (6) 

where  —<  u • u > /3  is  the  mean  square  of  one  component  of  the  velocity.  The 
proposal  is  to  numerically  solve 

^ + u- Vu  = -Vp/p  + W*u  + (e/3C/2)u  (7) 

with  the  objective  of  comparing  the  statistics  with  those  of  grid  flow  turbulence.  Equa- 
tion (7)  has  the  property  that  the  rest  state  is  unstable  to  long  waves.  Therefore  solutions 
cannot  decay  to  zero  but  must  transfer  energy  to  shorter  waves  in  order  to  dissipate  en- 
ergy. Some  reasons  are  given  below  for  expecting  computational  results  to  be  comparable 
with  experiments  for  inertial  scales  of  turbulence. 


1.2.  Linear  forcing  of  the  Karman-Howarth  equation 
A derivation  of  the  Karman-Howarth  equation,  following  the  steps  given  in  detail  by 
Landau  and  Lifshitz  (1959),  but  including  a forcing  term  as  in  (2),  gives 

1^_  1J__5  4 ^ 

dt  2 dt  6r^dr^  ^ 


^^<"(^2,«)-f(xi,t)  + u(xi,t)-f(x2,t)>dr  (8) 

where  B2(r,t)  and  Bs(r,t)  are  second  and  third  order  longitudinal  structure  functions, 

X2  = Xi  -j-  r and  r = |r|.  With  f given  by  (3),  the  last  term  can  be  written 

[ 2Qr^Rii{r,t)dr  (9) 

^ Jo 

where  Rij{r,t)  =<  Ui{xi,t)uj{xi  -f  r,<)  > is  the  velocity  correlation  tensor.  Since  the 
trace  of  the  correlation  tensor  is  related  to  the  second  order  structure  function  by 

= (10) 

(9)  can  be  integrated  to  the  form 

2QU^-QB2{r,t)  . (11) 


Using  this  with  Q given  by  (6),  and  dropping  the  time  derivatives  for  stationary  turbu- 
lence, gives  the  result 


e 

3U2 


J__5 

6r‘*  dr 


B2  = 


Id.  dB2 

u—T  — r 

dr  dr 


(12) 


in  which  e is  a constant  given  by  (5).  This  differs  from  the  standard  decaying  version  of 
the  Karman-Howarth  equation  which  is 


_2  _ 1^  ^ 1 4 _ ^dB2 

3^  2 dt  6r*dr^  ^ ^ r^  dr^  dr 

where  here  e is  a function  of  time  given  by 


(13) 


e(t)  - 


3 dU^  _ 15z/  (PB2 
2 dt  ~ 2 dr^ 


0 


(14) 
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The  Kolmogorov  law  may  be  derived  from  (12),  by  the  method  of  matched  asymp- 
totic expansions,  in  a manner  similar  to  that  employed  by  Lundgren  (2002)  to  get  that 
law  from  (13)  but  without  the  necessity  of  using  similarity  in  time.  The  result  is  the 
same: 


B2  = C^U^irJLf^^  = C2(er)2/3 

(15) 

where  L = U^/e  and  e and  are  independent  of  time  now. 

Equation(12)  can  be  integrated  to  obtain  Bz  in  terms  of  Bz,  as 

„ dBi  4 2e  1 r , 

(16) 

Using  (15)  for  Bz  gives  the  simple  result 

Bz  4C2/r\-i/3  6C2/r\5/3 

(17) 

This  may  be  rewritten  using  the  Taylor  microscale  A to  scale  r: 


Bs  = - Jer  ^1  - ^2.678C2(^)^^^  + 2.O29C2(0 

Here  Rl  = UL/v\  Ri  = -R^/15  relates  R\  and  Rl  and  X/L  = 15/J?a 
L.  The  corresponding  result  for  decaying  turbulence  (Lundgren  2002;  Lindborg  1999)  is 
almost  the  same: 

ft  = - fir'=(2.678(|  (») 

where  n is  the  energy  decay  exponent  oc  t~").  Equations  (18)  and  (19)  would  be 
exactly  the  same  if  n = 2,  which  is  not  a realistic  decay  exponent;  n = 1.2  is  often 
observed  and  n = 4/3  is  the  maximum  value  possible.  These  equations  give  a Reynolds 
number  correction  to  the  Kolmogorov  “4/5”  law  showing  that  it  is  approached  slowly  as 
Ra  —>■  00.  When  Cs  = 2 and  n = 1.2  the  compensated  forms  have  maxima  at  r/A  = 1.23 
(for  hnearly  forced  turbulence)  and  r/A  = 1.11  (for  decaying  turbulence). 

The  similarity  of  these  equations  can  be  seen  in  Figure  1 where  they  are  plotted  for 
= 350  and  compared  with  an  experiment  in  a turbulent  jet. 


j (18) 

relates  A and 


1.3.  Normalized  decaying  turbulence 

Consider  (7)  again,  with  and  without  the  isotropic  forcing  term.  In  the  stationary  case 
with  forcing  introduce  dimensionless  variables 


V = u/U,  X = r/L,  T = tU/L 

where  U and  L = /e  are  constants.  A simple  change  of  variables  gives 
|X+vVv  = -VP-l-i?7^V2u  + iv  . 

OT  3 

Now  consider  the  case  without  the  forcing  term,  with  the  change  of  variables 


X = r/Z.(t), 


dt 


(20) 


(21) 


V = u/U{t) 


(22) 
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Figure  1.  Compensated  third  order  structure  function  vs.  r/A  at  R\  = 350.  The  data  points  are 
from  a laboratory  jet  (Gagne,  2002)  at  the  same  value  of  R\.  The  long  dashed  curve  is  from  (19). 
The  dashed  curve  is  the  linear  forcing  result  from  (18).  The  solid  line  is  from  the  computation 
of  section  2 using  the  R^^^^  scaling  law  to  extrapolate  from  J?x  = 170  to  R>.  = 350. 

The  equation  transforms  to 

+ V • Vv  = -VP  + ~ • Vv.  (23) 

Using  £ = ^ — U^/e  the  coefficient  of  v is  so 

^+^r■'Vv  = -S/P  + Rl^V^^r+lv  + • Vv  . (24) 

C/i  t/ 

While  this  is  not  exactly  the  same  as  (21)  because  of  the  last  term,  it  is  apparent  that 
energy  decay  has  the  effect  of  an  isotropic  forcing  term.  Note  also  that  if  oc  t~V:  then 
e oc  and  L — /t  cx  So  L would  be  zero  if  n = 2.  In  this  case  the  equations 

would  be  the  same  except  for  the  time  dependence  of  P^- 

1.4.  Comparison  with  low  wavenumber  forcing 
In  this  subsection  a more  general  forcing  function  is  applied  to  the  Karman-Howarth 
equation  in  order  to  show  that  the  forcing  range  has  a significant  effect  on  the  third  or- 
der compensated  structure  function  and  presumably  on  any  inertial  range  statistics.  This 
analysis  is  carried  out  in  detail  in  Appendix  I and  only  summarized  here.  A Gaussian 
filter  is  applied  to  the  velocity  field,  filtering  out  a variable  part  of  the  high  wavenumber 
content.  This  filtered  velocity  is  used  as  a forcing  function  for  the  Karman-Howarth  equa- 
tion. In  one  limit  there  is  uniform  linear  forcing.  The  results  of  a calculation  are  shown 
in  Figure  2a  for  different  values  of  KL.  Here  K is  the  width  of  the  filter;  wavenumbers 
greater  than  K are  filtered  out.  When  K is  small  only  low  wavenumbers  are  forced,  while 
if  K — » 00  forcing  is  uniform  over  all  wavenumbers.  The  figure  shows  the  compensated 
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Figure  2a.  Third  order  compensated  structure  functions  for  variations  of  the  forcing  range  com- 
puted from  (A17).  The  curves  are  (from  the  bottom)  for  KL  = 1000, 100, 25, 10, 5, 3.  The  lowest 
curve  {KL  = 1000)  is  close  to  the  linear  forcing  result  from  (18).  All  the  curves  axe  computed 
for  Rx  = 1000. 


Figure  2b.  Solid  curves,  third  order  compensated  structure  function  vs.  r/rj  for  KL  = 5 and 
Rx  = 125,284,381,460  from  the  bottom  up.  This  compares  with  the  low  wavenumber  forced 
computation  of  Gotoh  et.  al.  (2002)  at  the  same  Reynolds  numbers  (compare  their  Fig.  12).  The 
maximum  values  shift  to  the  right  and  increase  with  increasing  Rx  in  a similar  manner.  The 
maxima  are  .70,  .76,  .77,  .78  comparing  with  .66,  .77,  .78,  .76  from  Gotoh.  The  dashed  curves 
are  for  linear  forcing  at  the  same  values  of  Rx  from  (18)  showing  the  differences. 
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Time 


Figure  3.  Time  History  of  R\  with  linear  forcing.  Numerical  data  for  the  following  figures  was 
taken  over  a time  range  of  about  .5  near  T=10,  where  R\  = 170. 


third  order  structure  function  for  different  values  oiKL  for  R\  = 1000.  When  KL  — 1000 
the  curve  is  nearly  the  same  as  computed  from  {18).  As  KL  decreases  the  maximum  value 
moves  to  the  right  and  increases  towards  4/5,  greatly  distorting  the  structure  function 
in  the  inertial  range. 

Figure  2b  shows  compensated  structure  functions  for  several  values  of  R\  for  the  single 
value  KL  = 5.  This  was  constructed  in  order  to  compare  with  the  high  resolution  DNS 
of  Gotoh  et.  al.  (2002),  which  was  driven  by  white  noise  in  a low  wavenumber  band 
which  corresponds  roughly  to  a Gaussian  filter  with  width  KL  = 5.  (In  their  numerics 
at  R\  = 460  L ~ 2.5,  making  K = 2 which  is  in  approximate  accord  with  their  forcing 
range  of  1 < A:  < \/6.) 


2.  Pseudospectral  Computation  with  Linear  Forcing 

Computations  have  been  carried  out  with  a Rogallo/Wray  pseudospectral  code,  mod- 
ified slightly  to  accomodate  linear  forcing.  This  was  done  with  (256)^  resolution  with 
u = .003,  box  size  27t  on  a side,  and  time  scale  set  by  taking  Q = 1.  Figure  3 shows 
the  time  history  of  Rx  during  a lengthy  run.  As  can  be  seen  the  computation  is  not 
exactly  stationary  (but  should  be  statistically  stationary  after  an  initial  transient).  The 
instantaneous  value  of  e/3f/^  (not  shown  here)  fluctuates  about  unity  by  about  ±20%. 
There  are  fairly  stationary  periods,  however.  Data  for  the  following  figures  was  taken  in 
a relatively  stationary  period  near  T = 10  where  R\  — 170. 

Figure  4 shows  the  energy  spectrum  versus  kt).  The  spectrum  is  a little  shallower  than 
^-5/3  upper  straight  line  on  the  figure).  Mydlarski  and  Warhaft  (1996)  indicate 
approximately  at  this  Reynolds  number  for  grid  flow  turbulence. 
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Figure  4.  Energy  spectrum  in  Kolmogorov  variables  for  Rx  — 170.  The  upper  straight  line  is 
and  the  lower  one  is 


Figure  5.  Compensated  even  order  structure  functions  versus  r/X  for  Rx  = 170.  Solid  lines  are 
curve  fit  using  (25).  The  solid  straight  lines  otc  an  extrapolation  to  Rx  = oo  using  V2,  V4,  Ve. 


In  Figure  5 compensated  structure  functions  of  second,  fourth,  and  sixth  order  are 
plotted  versus  r/X.  The  minimum  separation  between  two  points  (27t/256)  is  very  nearly 
.2A.  The  output  was  taken  by  averaging  over  all  pairs  of  points  along  the  x direction  with 
separation  r (which  is  a multiple  of  .2A).  For  each  structure  function  there  are  six  curves 
(the  lighter  long-dashed  curves  on  the  figure),  one  each  for  times  separated  by  1000  time 
steps.  The  total  elapsed  time  for  5000  time  steps  is  about  .5  time  units  on  Figure  3.  The 
spread  of  these  curves  indicates  the  random  nature  of  the  output  even  after  the  spatial 
averaging. 
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Figure  6.  Odd  order  compensated  structure  functions  for  R\  = 170.  Curve  fit  same  as  in  Figure 
5. 


The  pth  order  compensated  structure  function  should  take  the  form  (Lundgren  2003) 
5p/(er)P/3  = Fp  (^1  - + Bp(r/A)"V3)^  . (25) 

This  gives  Reynolds  number  corrections  to  Kolmogorov  (1941)  theory  and  can  be  re- 
garded as  a scaling  law.  If  the  parameters  Vp,  Ap,  Bp  were  determined  for  some  Reynolds 
number  (numerically  or  by  experiment)  one  could  extrapolate  the  data  to  another  Reynolds 
number.  The  coefficients  for  p = 2, 3, 4, 5, 6 were  determined  as  follows.  The  six  time  sets 
were  truncated  to  about  .5  < r/A  < 5 in  order  to  get  the  upper  parts  of  the  sets,  roughly 
centered  about  the  maxima.  These  were  accumulated  into  into  a single  set,  which  looks 
like  a lot  of  scattered  experimental  points.  Then  using  a nonlinear  curve  fitting  algorithm 
(on  xmgr)  the  coefficients  Vp,Ap,Bp  were  determined  for  R\  = 170  = .033).  The 

values  are  shown  in  the  table  below. 


p 

Fp 

Ap 

Bp 

(r/ A)max 

Vp/p\ 

2 

2.11 

2.62 

6.15 

2.17 

1.05 

3 

-.824 

6.56 

3.75 

1.07 

-.137 

4 

15.24 

4.85 

7.65 

1.78 

.635 

5 

-16.91 

7.81 

4.53 

1.08 

-.141 

6 

197.5 

6.56 

6.87 

1.45 

.274 

Table  1.  Coefficients  for  (25).  Fifth  column  is  position  of  maximum  of  Ra/er.  Sixth 
column  compares  Vp  with  p\. 


The  fit  curves  are  shown  as  heavy  lines  which  approximate  the  data  fairly  well.  These 
should  be  thought  of  as  time  averages  of  the  six  curves  of  each  set.  The  horizontal  straight 
lines  above  each  set  are  the  curves  5p/(er)P/^  = Vp,  representing  an  extrapolation  to 
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Figure  7.  Second  and  third  order  structure  functions  repeated  from  figures  5 and  6 on  a semilog 
plot  in  order  to  better  show  the  nature  of  the  curve  fit. 

R\  = oo.  In  particular  note  that  the  horizontal  line  above  the  second  order  stucture 
function,  V2  = 2.11,  is  an  acceptable  value  for  that  asymptote. 

The  third  and  fifth  order  structure  functions  were  treated  in  the  same  way  and  shown 
in  Figure  6.  The  individual  sets  show  much  more  scatter  than  the  even  order  structure 
functions  because  of  cancellation  between  the  left  and  right  sides  of  the  almost  symmetric 
pdfs.  Note  the  value  1^  = .824.  This  should  be  exactly  .8,  the  Kolmogorov  “4/5”  law. 
Figure  7 repeats  the  B2  and  83  curves  from  Figures  5 and  6 on  a semilog  plot. 

The  maxima  of  the  compensated  structure  functions  have  positions  which  scale  with  A. 
For  even  orders  the  positions  shift  towards  smaller  values  with  increasing  order.  For  odd 
orders  there  is  not  much  shift.  The  rapid  increase  of  the  magnitude  of  the  even  orders 
with  increasing  order  required  that  they  be  presented  on  a log-log  plot  in  order  to  show 
them  on  the  same  figure.  A rapid  increase  of  moments  like  pi,  as  seen  in  the  table,  is  a 
signature  for  exponential  tails  on  the  pdf.  This  is  observed  approximately  in  Figure  8. 

An  example  of  extrapolation  for  the  third  order  structure  function  is  shown  in  Figure 
1,  where  the  curve  fit  values  of  V3,  AsjBa  are  used  to  extrapolate  to  Rx  = 350  = 

.020).’ 

Figure  8 is  the  final  figure.  This  shows  the  velocity  pdf  versus  vja  where  v is  the  velocity 
difference  between  two  points  along  the  x direction  and  cr  is  the  standard  deviation, 
CT  z=<  In  these  variables  the  average  and  the  second  moment  are  both  unity. 

The  outermost  curve  is  as  close  to  the  pdf  of  the  velocity  derivative  as  can  be  reached 
with  this  resolution.  The  fluctuations  in  the  tails  were  ameliorated  to  some  extent  by 
accumulating  the  pdfs  over  10  successive  time  steps  (in  addition  to  spatial  averaging). 


3.  Conclusions 

There  were  two  objectives  to  this  research.  The  first  was  to  show  that  it  is  desirable  and 
useful  to  compute  stationary  isotropic  turbulence  with  a forcing  function  which  is  a simple 
linear  function  of  velocity,  which  thus  forces  uniformly  at  all  wavenumbers.  This  was  done 
by  analysis  of  of  a forced  Karman-Howarth  equation,  which  allows  calculation  of  third 
order  structure  functions  in  terms  of  second  order  structure  functions.  This  allowed  a 
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Figure  8.  pdf  of  velocity  difference  vs.  vfa,a  =<  >1/2  curves  are  for  values  of  the 

separation.  Prom  the  outside  rfX  = 2.,  A,  .6, 1.0, 1.4, 2.0, 2.4 


favorable  comparison  between  decaying  isotropic  turbulence  and  linearly  forced  isotropic 
turbulence  in  subsection  1.2. 

By  applying  filtered  linear  forcing  to  the  Karman-Howarth  equation  it  was  shown  that 
forcing  at  low  wavenumber  has  a large  effect  on  third-order  structure  functions,  and  would 
very  likely  influence  all  inertial  range  statistics,  in  partcular  it  could  affect  the  computa- 
tion of  anomalous  exponents.  The  third  order  structure  function  with  low  wavenumber 
filtered  forcing  was  compared  with  high  quality  DNS  (Gotoh  et.  al.,  2002)  and  linear 
forcing  was  compared  with  experiments  of  Gagne  (2002).  The  differences  between  low 
wave  number  forcing  and  linear  forcing  were  considerable,  as  seen  in  Figure  2b,  with  low 
wavenumber  forcing  approaching  the  four-fifths  law  too  rapidly  with  increasing  Reynolds 
number. 

The  second  objective  was  to  compute  a moderate  resolution  DNS  of  box  turbulence 
with  linear  forcing  to  show  the  feasibility  of  such  a computation.  This  was  done  in  section 
2.  Compensated  structure  functions  of  second  through  sixth  order  were  computed  at 
Rx  — 170.  It  was  shown  that  these  could  be  fit  with  an  scaling  law,  determining 

the  coefficients  to  fit  equation  (25).  (The  analysis  which  produced  this  law  (Lundgren 
2003)  applies  equally  to  the  forced  Navier-Stokes  equation.)  The  usefulness  of  the  scaling 
law  was  demonstrated  in  Figure  1,  where  the  third-order  structure  function  computed  at 
R\  — 170  was  extrapolated  to  R\  — 350  (the  solid  curve  in  Figure  8 was  extrapolated 
to  the  solid  curve  in  Figure  1),  where  it  compares  with  the  Karman-Howarth  results  and 
with  experiment. 
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5.  Appendix  I.  Filtered  Forcing. 

With  a general  forcing  term,  the  integrated  Karman-Howarth  equation  is 

B3  = - ^er  - / rdr  f < U2  • fi  + Ui  • £2  > r'^dr  . (Al) 

“r  5 r^  Jq  Jq 

It  is  desired  to  develop  a general  case  which  includes  forcing  at  small  wavenumbers 
and  also  broader  forcing.  A simple  model  takes  the  forcing  function  proportional  to  the 
velocity,  but  with  a variable  part  of  the  high  wavenumber  end  filtered  out  with  a Gaussian 
filter.  A filtered  velocity  in  Fourier  space  may  be  written 


u<(k)  = pjf(fc)u(k)  (A2) 

where  the  filter  and  its  transform  are 

gK{x)  = J exp(ik-x)gK(k)dk  (A3) 

^ ■ ^)9K{x)dx  (A4) 

It  is  assumed  that  (0)  = 1 so  that 

j 9K{x)dx  = (27t)^  . (A5) 

For  a Gaussian  filter 

gK  = exp{-.5k^/K^)  . (A6) 

This  filters  out  wavenumbers  /c  > A for  0 < A < 00;  A = 00  is  the  uniform  forcing  case. 
The  transform  of  gK  is 

p^(x)  = (277)3/2^3  exp(-AV/2)  . (A7) 

In  general  the  filtered  physical  space  velocity  is  given  by  a convolution: 

/ 5«-(|x-  sl)u(s)ds  . (A8) 

The  forcing  function  will  be  assumed  to  be  f = Qu'^,  where,  since  e =<  u • f > 


From  (A8) 


where 


Q = e/  <u-u^  > . 


< u 


>= 


(27t): 


/«(i 


X — s|)  < u(s)  • u(x)  > ds 


(A9) 

(AlO) 


< u(s)  • u(x)  >=  i2ii(|x  - s|)  (All) 

is  the  trace  of  the  correlation  function  given  by  (10).  By  using  (A9)  and  (AlO)  the  force 
correlation  in  (Al)  may  be  written 


< U2  • fi  + ui  • f2  >=  2e 


/gif(|x-s|)fiii(s)ds 
f gK{s)Rii{s)ds 


(A12) 


B3  — 


dB2 

dr 


Substituting  this  into  (Al)  gives 


(A13) 
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|er(27T)33t/^  + ^ rdr  r'^dr  j qk{\s  - r|)(  - \^!j;S^B2{s))ds 
{2n)W  + JgK{s){  - lJ,£s^B2{s))ds 

where  the  integration  over  the  3U^  part  of  Ru  was  done  using  (A5)  and  rdr  r'^dr  = 
r^/15.  In  the  upper  integral  do  the  integration  over  the  angle  variables  with  the  Gaussian 
filter,  obtaining 


/ 


where 


B{s,r)  = 


p/f(|s  — r\)  sm{&)ddd(p  — KB{s,r) 

f (exp(— .5AT^(s  — r)^)  — exp{—.5K^{s  + r)^)) 


sr 


(AU) 


(A15) 


is  a temporary  notation.  Now  do  the  integration  by  parts  on  s,  which  shifts  the  differen- 
tiation onto  B{s,r).  The  following  identity  may  be  proved 

dB{s,r)  _ 1 dA{s,r) 


ds 


K'^s'^r^  dr 


{A16) 


where 


A(s,  r)  = (1  — K'^sr)  exp(—.5K^{s  — r)^)  — (1  -1-  K'^sr)  exp(— .5K^(s  -f  r)^)  . 

Because  of  the  1/r^  factor  in  (A16)  the  inner  / r'^dr  may  be  carried  out,  resulting  in  the 
near-final  form 

Bz  _ 6L^  db2  I + (27r)-V^jr-^^  /J~  rdr  Ajs,  r)b2{s)sds 

er  Rl  dr  1 — {2'!t)~^/‘^K^j  exp{—.5K^s^)b2{s)s'^ds 

where  62  = ^2/1/^.  The  double  integrals  have  to  be  integrated  numerically  with  an 
assumed  function  for  B2.  An  appropriate  form  is 

B2  = 21/2  tanh(.5C'2(r/L)2/3)  (A18) 

which  gives  the  Kolmogorov  two-thirds  law  for  small  r/L  and  tends  to  the  proper  limit 
2172  as  r/L  — > 00. 
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